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Modelling Seismic Wave Propagation for 
Geophysical Imaging 


Jean Virieux et al.!* Vincent Etienne et al.2‘and Victor Cruz-Atienza et al.3# 
1ISTerre, Université Joseph Fourier, Grenoble 
2GeoAzur, Centre National de la Recherche Scientifique, Institut de Recherche 
pour le développement 
3 Instituto de Geofisica, Departamento de Sismologia, Universidad Nacional 
Autónoma de México 
12 France 
3 Mexico 


1. Introduction 


The Earth is an heterogeneous complex media from the mineral composition scale (~ 107 6m) 
to the global scale (~ 106m). The reconstruction of its structure is a quite challenging problem 
because sampling methodologies are mainly indirect as potential methods (Günther et al., 
2006; Rticker et al., 2006), diffusive methods (Cognon, 1971; Druskin & Knizhnerman, 1988; 
Goldman & Stover, 1983; Hohmann, 1988; Kuo & Cho, 1980; Oristaglio & Hohmann, 1984) or 
propagation methods (Alterman & Karal, 1968; Bolt & Smith, 1976; Dablain, 1986; Kelly et al., 
1976; Levander, 1988; Marfurt, 1984; Virieux, 1986). Seismic waves belong to the last category. 
We shall concentrate in this chapter on the forward problem which will be at the heart of any 
inverse problem for imaging the Earth. The forward problem is dedicated to the estimation 
of seismic wavefields when one knows the medium properties while the inverse problem is 
devoted to the estimation of medium properties from recorded seismic wavefields. 


The Earth is a translucid structure for seismic waves. As we mainly record seismic signals 
at the free surface, we need to consider effects of this free surface which may have a 
complex topography. High heterogeneities in the upper crust must be considered as well 
and essentially in the weathering layer zone which complicates dramatically the waveform 
and makes the focusing of the image more challenging. 


Among the main methods for the computation of seismic wavefields, we shall describe 
some of them which are able to estime the entire seismic signal considering different 
approximations as acoustic or elastic, isotropic or anisotropic, and attenuating effects. Because 
we are interested in seismic imaging, one has to consider methods which should be efficient 
especially for the many-sources problem as thousands of sources are required for imaging. 
These sources could be active sources as explosions or earthquakes. We assume that their 
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distribution are known spatially as punctual sources and that the source time function is the 
signal we need to reconstruct aside the medium properties. 


Asymptotic methods based on the high frequency ansatz (see (Virieux & Lambaré, 2007) 
for references or textbooks (Cerveny, 2001; Chapman, 2004)) and spectral methods based 
on a spatial and time Fourier transformations (Aki & Richards, 2002; Cagniard, 1962; de 
Hoop, 1960; Wheeler & Sternberg, 1968) are efficient methods which are difficult to control: 
whispering galeries for flat layers are efficiently considered using spectral methods. These 
two methods may be used either for local interpretation of specific phases or as efficient 
alternatives when media is expected to be simple. They could be used as well for scattering 
inverse problems. In the general heterogeneous case, we have to deal with volumetric 
methods where the medium properties are described through a volume while seismic wave 
fields satisfy locally partial differential equations. Although one may consider boundaries as 
the free surface or the water/solid interface, we may consider that variations of the medium 
properties are continuous at the scale of the wavelength which we want to reconstruct: 
the best resolution we could expect is half the wavelength (Williamson & Worthington, 
1993). Therefore a volumetric grid discretization is preferred where numerical expressions 
of boundary conditions should be mostly implicit through properties variations. 


A quite popular method because of this apparent simplicity is the finite difference method 
where partial derivatives are transformed into finite difference expressions as soon as the 
medium has been discretized into nodes: discrete equations should be exactly verified. We 
shall consider first this method as it is an excellent introduction to numerical methods and 
related specific features. We will consider both time and frequency approaches as they have 
quite different behaviours when considering seismic imaging strategies. 


Applications will enhance the different properties of this numerical tool and the caveats we 
must avoid for the various types of propagation we need. 


Another well-known approach is the finite element method where partial differential 
equations are asked to be fulfilled in a average way (to be defined) inside elements paving 
the entire medium. We shall concentrate into the discontinuous Galerkin method as it allows 
to mix acoustic and elastic wave propagation into a same formalism: this particular method 
shares many features of finite element formalism when describing an element, but differs by 
the way these elements interact each other. We avoid the description of the continuous finite 
element method for compactness and differences will be pointed out when necessary. Again, 
we shall discuss both time-domain and frequency-domain approaches. 


Applications will illustrate the different capabilities of this technique and we shall illustrate 
what are advantages and drawbacks compared to finite difference methods while specific 
features will be identified compared to continuous finite element methods. 


We shall conclude on the strategy for seismic imaging when comparing precision of solutions 
and numerical efforts for both volumetric methods. 


2. Equations of seismic wave propagation 


In a heterogeneous continuum medium, seismic waves verify partial differential equations 
locally. Integral equations may provide an alternative for the evolution of seismic fields either 
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in the entire domain or at the scale of an elementary element of a given mesh describing the 
medium structure. 


Fundamental laws of dynamics require the conservation of linear and angular momentum in 
a Galilean reference frame. In the continuum, a force applied across a surface, oriented by 
the unit normal n at a given position x = (x,y, z) in a Cartesian coordinate system (O, x,y, z), 
by one side of the material on the other side defines the traction vector t; = ojn; where the 
second-rank stress tensor 7 has been introduced. The conservation of the angular momentum 
makes the stress tensor symmetrical 0; = 0;;. We shall introduce as well a volumetric force 
per unit mass at the current position denoted as f = (fx, fy, fz). The conservation of linear 
momentum allows one to write the acceleration of the displacement motion u(x) of a given 
particle at the current position as 


du; dik 
Mes 


+ (x) fir (1) 


where the density is denoted by p(x). 


Aside the translation and the rotation transformations preserving the distances inside the 
body we consider, the deformation of the continuum body is described by defining a strain 


tensor € expressed as 
1 (uj Ou; 
,==(—4+ 2}. 2 
ij 2 (= + Ox; ( ) 


The symmetrical definition of the deformation ensures that no rigid-body rotations are 
included. The particle motion is decomposed into a translation, a rotation and a deformation: 
the two formers transformations preserve distances inside the solid body while the third one 
does not preserve distances, inducing stress variations inside the solid body. In the framework 
of linear elasticity, there is a general linear relation between the strain and stress tensors by 
introducing fourth-rank tensor c;;, defined as follows 


Tij = CijklE€kl- (3) 


Because of symmetry properties of stress and strain tensors, we have only 36 independent 
parameters among the 81 elastic coefficients while the positive strain energy leads to a further 
reduction to 21 independent parameters for a general anisotropic medium. For the particular 
case of isotropic media, we end up with two coefficients which can be the Lamé coefficients 
À and y. The second one is known also as the rigidity coefficient as it characterizes the 
mechanical shear mode of deformation. The following expression of elastic coefficients, 


Cijkl = Adjjôkr + (Sd; + 9:19jk), (4) 
with the Kronecker convention for 4;; gives the simplified expression linking the stress tensor 
to the deformation tensor for isotropic media as 


oj = NE KK ij + 2pE;;. (5) 
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One may prefer the inverse of the relation (5) where the deformation tensor is expressed from 
the stress tensor by the introduction of the Young modulus E. Still, we have two independent 
coefficients. By injecting the relation (5) into the fundamental relation of dynamics (1), we 
end up with the so-called elastic wave propagation system, which is an hyperbolic system of 
second order, where only the displacement u has to be found. This system can be written as 


Ou, 1 duy duy du, 
dE p a Pa geet tH) a T a | 


du 1 du du au 
y — | y | | x | Z | 
ae p a PA gg ATE (£ i S) | 
du du 
y y 
a | 0x2 0722 ) (6) 
du; 1 du, duy uy 
ae p a Page PAER a T ayaa | 


+ du, + du, 
PY ax dy? 7 
where we have neglected spatial variations of Lamé coefficients. Therefore, we must 
reconstruct over time the three components of the displacement or equivalently of the velocity 


or the acceleration. Choosing the stress is a matter of mechanical behaviour in a similar way 
for seismic instruments which record one of these fields. 


For heterogeneous media, spatial differential rules for Lamé coefficients have to be designed. 
We shall see how to avoid this definition in the continuum by first considering hyperbolic 
system of first-order equations, keeping stress field. More generally, any hyperbolic equation 
with n-order derivatives could be transformed in a hyperbolic system with only first 
derivatives by adding additional unknown fields. This mathematical transformation comes 
naturally for the elastodynamic case by selecting the velocity field v and the stress field 7 as 
fields we want to reconstruct. In a compact form, this first-order system in particle velocity 
and stresses is the following 


00; 
PJ = ei tel (7) 
a = Nog 5 0;; ji : 7b 
DE — Uk,kôij + u(vij F dji), (7b) 


with i,j = x,y,z. We may consider other dual quantities as (displacement, integrated stress) 
or (acceleration, stress rate) as long as the medium is at rest before the dynamic evolution. Let 
us underline that time partial derivatives are on the left-hand side and that spatial variations 
and derivations are on the right-and side. 
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Using simple linear algebra manipulations, alternative equivalent expressions may deserve 
investigation: the three components 0; could be linearly combined for three alternative 
components considering the trace o} = trace(æ)/3, the x-coordinate deviatoric stress o) = 
(20xx — Oyy — zz) /3 and the y-coordinate deviatoric stress 03 = (—Oxx + 20yy — zz) /3 which 
allows to separate partial spatial derivatives in the right hand side and material properties in 
the left hand side. The system (7) becomes 


dv; 
Pa = Mj + Of 
+ i 
BA+2n ot 
3 002 _ [00x 7 
2u at (3 ax vi) (8) 
3 003 o 00y 
2u at (3 oy ae 
1 00; 
mot it ji 


which could be useful when we move from differential formulation to integral formulation 
over elementary volumes. Partial differential operators only in the right-hand side of the 
system (8) are separated from spatial variations of model parameters on the left-hand side 


as a diagonal matrix À = ( sy oy ay 7 7 3) Similar strategies could be applied for 2D 
geometries. 


Finally, for easing discussions on the numerical implementation, let us write both the 1D scalar 
second-order acoustic wave equation in the time domain as 


u(x,t) 9 du(x, t) 


p(x) at2 = JEU) ax 7 (9) 
or, in frequency domain, 
0 du(x,w) 
2 | z = 
w o(x)u(x,w) + JEU) R 0, (10) 


away from sources where one can see the importance of the mixed operator 0,E(x)dx. 
We have introduced the Young modulus E related to unidirectional compression/delation 
motion. The 1D vectorial first-order acoustic wave equation can be written as 


mes 200 _ 2 t) 


d(x, t) dv(x, t) 
at OE 


(11) 
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from which one can deduce immediatly the system of equations in the frequency domain 


—iwp(x)o(x,w) = 


dv(x,w) 


—iwo(x,w) = E(x) 3x 


(12) 
Please note that the mixed operator does not appear explicitly. By discretizing this system and 
by eliminating the stress discrete values, one can go back to an equation involving only the 
velocity: a natural and systematic procedure for discretizing the mixed operator as proposed 
by Luo & Schuster (1990). 


For an isotropic medium, two types of waves - compressional and shear waves - are 
propagating at two different velocities vp and vs. These velocities can be expressed as 


Up = Ata and m = jë, (13) 


except for the 1D medium where only compression/dilatation motion could take place. The 
displacement induced by these two different waves is such that compressive waves uP verify 
V x uP = 0 and shear waves u° verify V - u$ = 0. Applying these operators to the numerical 
displacement will separate it into these two wavefields. 


2.1 Time-domain or frequency-domain approaches 


These systems of equations could be solved numerically in the time domain or in the 
frequency domain depending on applications. For seismic imaging, the forward problem 
has to be solved for each source and at each iteration of the optimisation problem. The 
time approach has a computational complexity increasing linearly with the number of 
sources while precomputation could be achieved in the frequency domain before modelling 
the propagation of each source. Let us write a compact form in order to emphasize the 
time/frequency domains approaches. The elastodynamic equations are expressed as the 
following system of second-order hyperbolic equations, 


w(x, t) 


MG)— 5 


= S(x)w(x,t) + s(x,t), (14) 
where M and S are the mass and the stiffness matrices (Marfurt, 1984). The source term is 
denoted by s and the seismic wavefield by w. In the acoustic approximation, w generally 
represents pressure, while in the elastic case, w generally represents horizontal and vertical 
particle displacements. The time is denoted by t and the spatial coordinates by x. Equation 
(14) is generally solved with an explicit time-marching algorithm: the value of the wavefield 
at a time step (n + 1) at a spatial position x is inferred from the value of the wavefields at 
previous time steps (Dablain, 1986; Tal-Ezer et al., 1990). Implicit time-marching algorithms 
are avoided as they require solving a linear system (Marfurt, 1984; Mufti, 1985). If both 
velocity and stress wavefields are helpful, the system of second-order equations can be recast 
as a first-order hyperbolic velocity-stress system by incorporating the necessary auxiliary 
variables (Virieux, 1986). The time-marching approach could gain in efficiency if one consider 
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local time steps related to the coarsity of the spatial grid (Titarev & Toro, 2002): this leads 
to a quite challenging load balancing program between processors when doing parallel 
programming as most processors are waiting for the one which is doing the maximum of 
number crunching as illustrated for the ADER scheme (Dumbser & Käser, 2006). Adapting the 
distribution of the number of nodes to each processor depending on the expected complexity 
of mathematical operations is still an open problem. Other integration schemes as the 
Runge-Kutta scheme or the Stormer/Verlet symplectic scheme (Hairer et al., 2002) could be 
used as well. 


Seismic imaging requires the cross-correlation in time domain or the product in frequency 
domain of the incidents field of one source and the backpropagated residues from the receivers 
for this source. In order to do so, one has to save at each point of the medium the incident 
field from the source which could be a time series or one complex number. The storage when 
considering a time-domain approach could be an issue: a possible strategy is storing only 
few time snapshots for recomputing the incident field on the fly (Symes, 2007) at intermediate 
times. An additional advantage is that the attenuation effect could be introduced as well. in 
the time-domain approach, the complexity increases linearly with the number of sources. 


In the frequency domain, the wave equation reduces to a system of linear equations, the 
right-hand side of which is the source, and the solution of which is the seismic wavefield. 
This system can be written compactly as 


B(x,w)w(x,w) = s(x,w), (15) 


where B is the so-called impedance matrix (Marfurt, 1984). The sparse complex-valued 
matrix B has a symmetric pattern, although is not symmetric because of absorbing boundary 
conditions (Hustedt et al., 2004; Operto et al., 2007). The fourier transform is defined with the 
following convention 


flo) = [fetta 


Solving the system of equations (15) can be performed through a decomposition of the matrix 
B, such as lower and upper (LU) triangular decomposition, which leads to the so-called 
direct-solver techniques. The advantage of the direct-solver approach is that, once the 
decomposition is performed, equation (15) is efficiently solved for multiple sources using 
forward and backward substitutions (Marfurt, 1984). This approach has been shown to be 
efficient for 2D forward problems (Hustedt et al., 2004; Jo et al., 1996; Stekl & Pratt, 1998). 
However, the time and memory complexities of the LU factorization, and its limited scalability 
on large-scale distributed memory platforms, prevents the use of the direct-solver approach 
for large-scale 3D problems (i.e., problems involving more than ten millions of unknowns) 
(Operto et al., 2007). 


Iterative solvers provide an alternative approach for solving the time-harmonic wave equation 
(Erlangga & Herrmann, 2008; Plessix, 2007; Riyanti et al., 2006; 2007). Iterative solvers are 
currently implemented with Krylov-subspace methods (Saad, 2003) that are preconditioned 
by the solution of the dampened time-harmonic wave equation. The solution of the dampened 
wave equation is computed with one cycle of a multigrid. The main advantage of the iterative 
approach is the low memory requirement, while the main drawback results from the difficulty 
to design an efficient preconditioner, because the impedance matrix is indefinite. To our 
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knowledge, the extension to elastic wave equations still needs to be investigated. As for the 
time-domain approach, the time complexity of the iterative approach increases linearly with 
the number of sources or, equivalently, of right-hand sides, in contrast to the direct-solver 
approach. 


An intermediate approach between the direct and the iterative methods consists of a hybrid 
direct-iterative approach that is based on a domain decomposition method and the Schur 
complement system (Saad, 2003; Sourbier et al., 2011): the iterative solver is used to solve the 
reduced Schur complement system, the solution of which is the wavefield at interface nodes 
between subdomains. The direct solver is used to factorize local impedance matrices that are 
assembled on each subdomain. Briefly, the hybrid approach provides a compromise in terms 
of memory saving and multi-source-simulation efficiency between the direct and the iterative 
approaches. 


The last possible approach to compute monochromatic wavefields is to perform the modeling 
in the time domain and extract the frequency-domain solution, either by discrete Fourier 
transform in the loop over the time steps (Sirgue et al., 2008) or by phase-sensitivity detection 
once the steady-state regime has been reached (Nihei & Li, 2007). An arbitrary number of 
frequencies can be extracted within the loop over time steps at a minimal extra cost. Time 
windowing can easily be applied, which is not the case when the modeling is performed in the 
frequency domain. Time windowing allows the extraction of specific arrivals (early arrivals, 
reflections, PS converted waves) for the full waveform inversion (FWI), which is often useful 
to mitigate the nonlinearity of the inversion by judicious data preconditioning (Brossier et al., 
2009; Sears et al., 2008). 


Among all of these possible approaches, the iterative-solver approach has theoretically the 
best time complexity (here, complexity denotes how the computational cost of an algorithm 
grows with the size of the computational domain) if the number of iterations is independent 
of the frequency (Erlangga & Herrmann, 2008). In practice, the number of iterations generally 
increases linearly with frequency. In this case, the time complexity of the time-domain 
approach and the iterative-solver approach are equivalent (Plessix, 2007). 


For one-frequency modeling, the reader is referred to those articles (Plessix, 2007; 2009; 
Virieux et al., 2009) for more detailed complexity analysis of seismic modeling based on 
different numerical approaches. A discussion on the pros and cons of time-domain versus 
frequency-domain seismic modeling relating to what it is required for full waveform inversion 
is also provided in Vigh & Starr (2008) and Warner et al. (2008). 


2.2 Boundary conditions 


In seismic exploration, two boundary conditions are implemented for wave modeling: 
absorbing boundary conditions to mimic an infinite medium and free surface conditions on 
the top side of the computational domain to represent the air-solid or air-water interfaces 
which have the highest impedance contrast. For internal boundaries, we assume that effects 
are well described by variations of the physical properties of the medium: the so-called 
implicit formulation (Kelly et al., 1976; Kummer & Behle, 1982). 
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2.2.1 PML absorbing boundary conditions 


For simulations in an infinite medium, an absorbing boundary condition needs to be applied 
at the edges of the numerical model. An efficient way to mimic such an infinite medium 
can be achieved with Perfectly-Matched Layers (PML), which has been initially developed 
by Berenger (1994) for electromagnetics, and adapted for elastodynamics by Chew & Liu 
(1996); Festa & Vilotte (2005). PMLs are anisotropic absorbing layers that are added at the 
periphery of the numerical model. The classical PML formulation is based on splitting of 
the elastodynamic equations. A new kind of PML, known as CPML, does not require split 
terms. The CPML originated from Roden & Gedney (2000) for electromagnetics was applied 
by Komatitsch & Martin (2007) and Drossaert & Giannopoulos (2007) to the elastodynamic 
system. CPML is based on an idea of Kuzuoglu & Mittra (1996), who has obtained a strictly 
causal form of PML by adding some parameters in the standard damping function of Berenger 
(1994), which enhanced the absorption of waves arriving at the boundaries of the model with 
grazing incidence angles. 


In the frequency domain, the implementation of PMLs consists of expressing the wave 
equation in a new system of complex-valued coordinates ¥ defined by (e.g., Chew & Weedon, 
1994): 
0 1. O 
— Ez) ax (16) 
In the PML layers, the damped 1D acoustic wave equation could be deduced from the 
equation (10) as 
1. 0 E(x) à 
2 j 
&w°p(x) + 
p(x) Ex(x) Ax x(x) ax 
where y(x) = 1+ i7x(x)/w and x(x) is a 1D damping function which defines the PML 
damping behavior in the PML layers. In the CPML layers, the damping function ¢,(x) 
becomes 


u(x,w) = —s(x,w), (17) 


dx 
, = 1 
Gx(x) = tx + ay + iw” (18) 


with angular frequency w and coefficients x; > 1 and «>; > 0. The damping profile dy varies 
from 0 at the entrance of the layer, up to a maximum real value demay at the end (Collino & 
Tsogka, 2001) such that 


by 4? 
dx = denas ( g5) : (19) 
and 
log (Reveff) 
Axmax = —3Vp OL coeff 2 (20) 
cpml 


with ôx as the depth of the element barycentre inside the CPML, Lepml the thickness of the 
absorbing layer, and Reoeff the theoretical reflection coefficient. Suitable expressions for 
Kx, dy and a are discussed in Collino & Monk (1998); Collino & Tsogka (2001); Drossaert 
& Giannopoulos (2007); Komatitsch & Martin (2007); Kuzuoglu & Mittra (1996); Roden & 
Gedney (2000). We often choose Rege¢¢ = 0.1% and the variation of the coefficient a; goes 
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from a maximum value (4,,,, = Æfo) at the entrance of the CPML, to zero at its end. If xy = 1 
and a, = 0, the classical PML formulation is obtained. 


One can use directly these frequency-dependent expressions when considering the frequency 
approach. The formulation in the time domain is slightly more involved. The spatial 
derivatives are replaced by 


1 
Ox => 0x + Ga * Ox, cn 
x 
with 
x(t) = -E ppe tant, i 
x 


where H(t) denotes the Heaviside distribution. Roden & Gedney (2000) have demonstrated 
that the time convolution in equation (21) can be performed in a recursive way using memory 
variables defined by 


Px = Gx * Ox. (23) 


The function Yx represents a memory variable in the sense that it is updated at each time 
step. Komatitsch & Martin (2007) have shown that the term xy has a negligible effect on the 
absorbing abilities, and it can be set to 1. If we take x, = 1, we derive the equation (23) using 
the equation (22) as 


Ory = —dydy — (dx + Ay) Wy. (24) 


One equation is generated for each spatial derivative involved in the elastodynamic system, 
which can be a memory-demanding task. Once they are computed at each time step, we 
can introduce the memory variables into the initial elastodynamic system which requires two 
additional variables for the 1D equations (11) with the definition of ~,(v) and ,(c) leading 
to the following system 


PDT = E + Vale) 
be = EE + px(0) ey 
to) = Toen (d(x) + ax(x))p(v) 
dpx(c) = op(c) | 
D = -dy (x)= — Gel) + az(x))p(o) 


At the outer edge of the PML zone, one could apply any conditions as simple absorbing 
conditions (Clayton & Engquist, 1977) or free surface conditions (Etienne et al., 2010) as fields 
go to zero nearby the outer edge. 


We must underline that the extension to 2D and 3D geometries is straightforward both in the 
frequency domain (Brossier et al., 2010; 2008) and in the time domain (Etienne et al., 2010; 
Komatitsch & Martin, 2007). 
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2.2.2 Free surface 


Planar free surface boundary conditions can be simply implemented using a strong 
formulation or a weak formulation. 


In the first case which is often met in the finite-difference methods, one requires that the stress 
is zero at the free surface. The free surface matches the top side of the FD grid and the stress 
is forced to zero on the free surface (Gottschamer & Olsen, 2001). Alternatively, the method 
of image can be used to implement the free surface along a virtual plane located half a grid 
interval above the topside of the FD grid (Virieux, 1986). The stress is forced to vanish at the 
free surface by using a virtual plane located half a grid interval above the free surface where 
the stress is forced to have opposite values to that located just below the free surface. In case 
of more complex topographies, one strategy is to adapt the topography to the grid structure 
at the expense of numerical dispersion effect (Robertsson, 1996) or to deform the underlying 
meshing used in the numerical method to the topography (Hestholm, 1999; Hestholm & Ruud, 
1998; Tessmer et al., 1992). In the first case, because of stair-case approximation, a local fine 
sampling is required (Hayashi et al., 2001). 


Owing to the weak formulation used in finite-element methods, the free surface boundary 
condition are naturally implemented by zeroing test functions on these boundaries which 
follow edges of grid elements (Zienkiewicz & Taylor, 1967). This approach could be used as 
well for finite-difference methods through the summation-by-parts (SBP) operators based on 
energy minimization combined with Simultaneous Approximation Term (SAT) formulation 
based ona boundary value penalty method (Taflove & Hagness, 2000). In this case, boundaries 
on which stress should be zero are not requested to follow any grid discretisation. 


Finally, one may consider an immersed boundary approach where the free surface boundary 
is not related to the discretisation of the medium as promoted by LeVeque (1986). Extensive 
applications have been proposed by Lombard & Piraux (2004); Lombard et al. (2008) where 
grid discretisation does not influence the application of boundary conditions. This approach 
might be seen as an extension of the method of images following extrapolation techniques 
above the free surface at any a priori order of precision. 


2.3 Source implementation 


There are different ways of exciting the numerical grid by the source term. The simplest one 
is the direct contribution of the source term in the discrete partial differential equations: for 
example, we may just increment by the source term after each time step or we may consider 
the right-hand side source term for solving the linear system in the frequency domain. 
Depending on the numerical approach, it is necessary to consider specific influences coming 
from the discretisation as we shall see for numerical methods we consider. 


In order to avoid singularities of solutions nearby the source, on can use the injection 
technique as proposed in the pioneering work of Alterman & Karal (1968) where a specific 
box around the source is defined. Inside the box, only the scattering field is computed. The 
incident field is estimated at the edges of the box and it is substracted when propagations 
are estimated inside the box and added when propagations are estimated outside the box. 
A more general framework is proposed by Oprèal et al. (2009) related to boundary integral 
approaches. 


264 Seismic Waves, Research and Analysis 


3. Finite-difference methods: solving the equations through a strong formulation 


We shall first consider the discretization based on simple and intuitive approaches as 
finite-difference methods for solving these partial differential equations while focusing on 
techniques useful for seismic imaging which means a significant number of forward problems 
for many sources in the same medium. We shall identify features which might be interesting 
for seismic imaging. If these approaches are intuitive for solving differential equations, the 
numerical implementation of boundary conditions and source excitation is less obvious and 
requires specific strategies as we shall see. 


3.1 Spatial-domain finite-difference approximations 


Whatever is our strategy for the reconstruction of the wave field u, one has to discretize it. 
We may be very satisfied by considering a set of discrete values (u4, U2,...,Uj_1,Uy) along one 
direction at a given specific time tn which can be discretized as well. Therefore, a simple way 
of solving this first-order differential system is by making finite difference approximations of 
spatial derivatives. 


Still considering a 1D geometry, the partial operator (0/dx) could be deduced from a Taylor 
expansion using Lagrange polynomial. A quite fashionable symmetrical estimation using a 
centered finite difference approximation is expressed as 

duf Uy — Up 


T 2Ax +0 [a+] ; 6) 


which is a three-nodes stencil as three nodes are involved: two for the derivative estimation 
and one for the updating. Let us mention that the discrete derivative is shifted with respect 
to discrete values of the field. Because of the very specific antisymmetrical structure of our 
first-order hyperbolic system where time evolution of velocity requires only stress derivatives 
(and vice versa), we may consider centered approximations both in space and in time. This 
will lead to a leap-frog structure or a red/black pattern. Of course, we have truncation errors 
expressed by the function O(Ax”) which depends on the power n of the spatial stepping and 
by the function O(At*) on the power k of the time stepping. 


We may require a greater precision of the derivative operator by using more points for this 
partial derivative approximation and a very popular centered finite difference approximation 
of the first derivative is the following expression 


n n n n 
ðu” “1 (aa T ut) Te Ce g ut) 
i — 2 2 2 2 j 4 
ox Ax | Ofax | a 
where cy = 9/8 et c2 = —1/24 (Levander, 1988). This fourth-order stencil is compact 


enough (few discrete points inside the stencil) for numerical efficiency while having a small 
local truncation error. This stencil is a five-nodes stencil. Let us underline that centered 
approximations lead to have field quantities not at the same position in the numerical grid as 
derivative approximations (figure 1). In other words, stress and velocity components should 
be specified on different positions of the spatial grid. If we still consider a full grid where stress 
and velocities are known at the same position, this stencil could be recast as a seven-nodes 
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Fig. 1. Cell structure for the full staggered grid (top), the partial staggered grid (bottom left) 
and the full grid (bottom right, Stress tensor is denoted by øg, particle velocity by v and the 
density by p as well as Lamé coefficient by À et u. The regular grid step is denoted by h. Time 
stepping is missing 


stencil given by 


n n n n n n 
ðu! di (ut, z ut) + do (ut. =- ul! 5) + ds (ut, =- u 3) 


1 


ox Ax! g 


(28) 


where Ax’ = Ax/2 and where we have following specific coefficients dj = c1, dy = 0. 
and d3 = c2. The fourth-order scheme would require the following theoretical coefficients 
dı = 15/20, dy = —3/20 and d = 1/60. For fourth-order stencils, the two sub-grids 
are not entirely decoupled and are weakly coupled leading to a dispersion behaviour as 
if the discretization is Ax. Let us remind that these sub-grids are completely decoupled 
when considering second-order stencils, leading to the staggered structure. Therefore, solving 
partial differential equations in the staggered grid structure has a less accurate resolution but 
improves significantly the efficiency of the method than solving equations in the full grid even 
with dispersion-relation-preserving stencils (Tam & Webb, 1993). The memory saving can be 
easily seen when comparing nodes for staggered grid and nodes for full grid (figure 1) 


When dealing with 2D and 3D geometries, we may exploit the extra freedom and estimate 
derivatives along the direction x from nodes shifted by half the grid step in x but also by 
half the grid step in y (and eventually in z). This leads to another compact stencil as shown 
in the figure 1 where all components of the velocity are discretized in one location while all 
components of the stress field are discretized half the diagonal of the grid as proposed by 
Saenger et al. (2000). This grid is still partially staggered and could be named as a partial grid. 


These standard and partial staggered structures are sub-grids of the full grid as shown in the 
figure 1 which is used in aeroacoustics (Tam & Webb, 1993). 
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These different discretizations related to various stencils may lead to preferential directions 
of propagation. Numerical anisotropy effect is observed even when considering isotropic 
wave propagation. The figure 2 shows error variations in velocities with respect to angles 
of propagation for the standard grid and the partial grid: one can see that the anisotropy 
behavior is completely different with a rotation shift of 45°. In 2D, the two stencils provide 
the same anisotropic error while the partial grid has a slightly improved numerical anisotropic 
performance (percentage differences go from 3 % down to 2 % in 3D geometries). Of course, 
the spatial sampling is such that the error should be negligible and few percentages is 
considered to be acceptable except nearby the source. 


Y 


Fig. 2. Numerical anisotropic errors when considering finite difference stencils related to 
partial staggered grid (left panel) and standard staggered grid (right panel) Saenger et al. (2000). 


Other spatial interpolations are possible. Previous discrete expressions are based on Lagrange 
interpolations while other interpolations are possible such as Chebychev or Laguerre 
polynomial or Fourier interpolations (Kosloff & Baysal, 1982; Kosloff et al., 1990; Mikhailenko 
et al., 2003). Interpolation basis could be local (Lagrange) or global(Fourier) ones based 
on equally spaced nodes or judiciously distributed nodes for keeping interpolation errors 
as small as possible: this will have a dramatic impact on the accuracy of the numerical 
estimation of the derivative and, therefore, on the resolution of partial differential equations. 
We should stress that local stencils should be preferred for seismic imaging for efficiency in 
the computation of the forward modeling. 


3.2 Time-domain finite-difference approximations 


Similarly, one may consider finite difference approximation for time derivatives which can be 
illustrated on the simple scalar wave equation. A widely used strategy is again the centered 
differences through the expression 


ðu” yori _ yt) 
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For understanding how the procedure of computing new values in time is performing, let us 
consider the simple 1D second-order scalar wave equation for displacement u. This equation 


u Pu 


I age (80) 
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could be discretized through these finite difference approximations 


41 = 
u Qup tu | 2 js — 2u; + “a (31) 


(At)? | (Ax)? 


The next value at the discrete time n + 1 comes from older values known at time n and time 
n — 1 through the expression 


u”, — 2u” + u” 
+1 wy 2 i+1 i i—1 fi —1 
ur & (cAt) | (Ax)? | tue ur. (32) 
A more compact notation of this equation as 
ee = 2(1- Shut + Se eae + ui) — yet (33) 
shows the quantity 
_ CA 
7 Ax’ 


known as the Courant number in the literature. This quantity is quite important for 
understanding the numerical dispersion and stability of finite difference schemes. The related 
stencil on the spatio-temporal grid as shown in the left panel of the figure 3 clearly illustrates 
that the value at time n + 1 is explicitly computed from values at time n — 1 and time n. In this 
explicit formulation, the selection of the time step At should verify that the Courant number 
is lower than 1 for any point of the medium. 


n+1 


i-1 i i+1 x i-1 i i+] x 


Fig. 3. Space/time finite difference stencil for an explicit scheme on the left and for an 
implicit scheme in a 1D configuration: black circles are known values from which the white 
circle is estimated. 


On the contrary, we may consider spatial derivatives at time n + 1. This leads us to an implicit 
scheme where more than one value at time n + 1 is present in the discretisation. The equation 
is now 

gett — 2u} + wet 8 Wu = au + yeu 


(At)? i (Ax)? 


(34) 
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which can be described by the Courant number S through the equation 


(1+ aS art? — S (uri + utt) = 2u} — ut, (85) 
The right panel of the figure 3 illustrates the structure of the stencil and that three unknowns 
have to be estimated through the single equation (35). By considering different spatial nodes, 
we may find these three unknowns by solving a linear system. The Courant number could 
take any value for time integration as long as discrete sampling is correctly performed. 


Other implicit stencils might be designed by averaging the spatial derivatives over the three 
times n — 1, n and n + 1. We may as well average the time derivative over the three positions 
i—1,iand i+ 1. This lead to another linear system to be solved. These weighting strategies 
could be designed for reducing numerical noise as numerical dispersion and/or anisotropy: 
a road for further improvements. 


As discretisation in space and time goes to zero, one expects the solution to be more precise 
but cumulative rounding errors should prevent to have too small values. In expressions (26) 
and (29), truncation error O[Ax?] goes to zero as the square of the discrete increment. We 
shall say that this is a second-order precision scheme both in space and in time. One consider 
often the stencil with the fourth-order precision in space and second-order precision in time, 


denoted as O| Ax*, At? |, as an optimal one for finite-differences simulations. 
p 


3.3 Frequency-domain finite-difference approximations 


The second-order acoustic equation (10) provides a generalization of the Helmholtz equation. 
In exploration seismology, the source is generally a local point source corresponding to an 
explosion or a vertical force. 


Attenuation effects of arbitrary complexity can be easily implemented in equations (10) and 
(12) using complex-valued wave speeds in the expression of the bulk modulus, thanks to 
the correspondence principle transforming time convolution into products in the frequency 
domain: in the frequency domain, one has to replace elastic coefficients by corresponding 
viscoelastic complex moduli for considering visco-elastic behaviors (Bland, 1960). For 
example, according to the Kolsky-Futterman model (Futterman, 1962; Kolsky, 1956), the 
complex wave speed ¢ is given by 


-1 
é=c | (2+ loseer +e) , (36) 


where the P wave speed is denoted by c = y E/p, the attenuation factor by Q and a reference 
frequency by w,. The function sgn gives the sign of the function. 


Since the relationship between the wavefields and the source terms is linear in the first-order 
and second-order wave equations, one can explicitely expressed the matrix structure of 
equations (10) and (12) through the compact expression, 


[M+S]u=Bu=s, (37) 
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where M is the mass matrix, S is the complex stiffness/damping matrix. This expression 
holds as well in 2D and 3D geometries. The dimension of the square matrix B is the number 
of nodes in the computational domain multiplied by the number of wavefield components. 
System (37) could be solved using a sparse direct solver. A direct solver performs first a LU 
decomposition of B followed by forward and backward substitutions for the solutions (Duff 
et al., 1986) as shown by the following equations: 


Bu = (LU)u=s (38) 
Ly=s; Uu=y (39) 


Exploration seismology requires to perform seismic modeling for a large number of sources, 
typically, up to few thousands for 3D acquisition. The use of direct solver is the efficient 
computation of the solutions of the system (37) for multiple sources. Combining different 
stencils for constructing a compact and accurate stencil can follow strategies developped for 
acoustic and elastic wave propagation (Jo et al., 1996; Operto et al., 2007; Stekl & Pratt, 1998). 
The numerical anisotropy is dramatically reduced 


The mass matrix M is a diagonal matrix although never explicitly constructed when 
considering explicit time integration. In the frequency domain formulation, we may spread 
out the distribution of mass matrix over neighboring nodes in order to increase the precision 
without increasing the computer cost as we have to solve a linear system in all cases. This 
strategy is opposite to the finite element approach where often the mass matrix is lumped into 
a diagonal matrix for explicit time integration (Marfurt, 1984). For a frequency formulation, 
considering the mass matrix as a non-diagonal matrix does not harm the solver. The weights 
of distribution are obtained through minimization of the phase velocity dispersion in an 
infinite homogeneous medium (Brossier et al., 2010; Jo et al., 1996): the numerical dispersion 
is dramatically reduced. 


3.4 PML absorbing boundary condition implementation 


Implementation of PML conditions in the frequency domain is straightforward using 
unsplit variables while, in the time domain, we need to introduce additional variables for 
handling the convolution through memory variables or to use split unphysical field variables 
(Cruz-Atienza, 2006). These additional variables are only necessary in the boundary layers 
following the figure 4 


We first consider an infinite homogeneous medium which is embedded into a cubic box of a 
16 km size and a grid stepping of h = 100 m. The thickness of the PML layer is 1 km leading 
to nsp = ten nodes inside the PML zone. The P-wave velocity is 4000 m/s while the S-wave 
velocity is 2300 m/s and the density 2500 kg/m°. The figure 5 shows various time sections of 
the 3D volume for the vertical particle velocity where one can see that the explosive wavefront 
is entering the PML zone at the time 2.8 s. The last two snapshots shows the vanishing of the 
wavefront with completely negligible residues at the final time (the decrease of the elastic 
energy is better than 0.2 % for ten nodes and could reach 0.03 % for twenty nodes). 


When we have discontinuous interfaces crossing the PML zone, we may expect difficulties 
coming from various angles of propagation waves (Chew & Liu, 1996; Festa & Nielsen, 2003; 
Marcinkovich & Olsen, 2003). Therefore, a simple heterogeneous medium is considered 
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nsp z 


Fig. 4. Three kinds of PML boundary layers should be considered where only one coordinate 
is involved (yellow zone), two coordinates are involved (blue zone) and three coordinates are 
involved (red zone). Internally, standard elastodynamic equations are solved 
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Fig. 5. Snapshots for y=0 of the vertical particle velocity at different times for an explosive 
source: on the left for an homogeneous infinite medium and on the right for an 
heterogeneous medium. Please note the vanishing of the seismic waves, thanks to the PML 


absorption 
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with two layers where physical parameters are (Up, Vs, 0)=(4330 m/s,2500 m/s, 2156 kg/ m?) 
and (Vp, Vs, p)=(6000 m/s,4330 m/s,2690 kg/ m3). The figure 5 shows that, in spite of the 
complexities of waves generated at the horizontal flat interface, the PML layer succeeds to 
absorb seismic energy with a residual energy of 0.3 % in this case, still far better than standard 
paraxial absorbing boundary conditions (Clayton & Engquist, 1977). 


3.5 Source and receiver implementation on coarse grids 


Seismic imaging by full waveform inversion is initiated at an initial frequency as small as 
possible to mitigate the non linearity of the inverse problem resulting from the use of local 
optimization approach such as gradient methods. The starting frequency for modeling in 
exploration seismics can be as small as 2 Hz which can lead to grid intervals as large as 200 m. 
In this framework, accurate implementation of point source at arbitrary position in a coarse 
grid is critical. One method has been proposed by Hicks (2002) where the point source is 
approximated by a windowed Sinc function. The Sinc function is defined by 
sin(7x) 


sinc(x) = pa « (40) 


where x = (xg — Xs), Xg denotes the position of the grid nodes and xs denotes the position 
of the source. The Sinc function is tapered with a Kaiser function to limit its spatial support 
(Hicks, 2002) . For multidimensional simulations, the interpolation function is built by tensor 
product construction of 1D windowed Sinc functions. If the source positions matches the 
position of one grid node, the Sinc function reduces to a Dirac function at the source position 
and no approximation is used for the source positioning. If the spatial support of the Sinc 
function intersects a free surface, part of the Sinc function located above the free surface is 
mirrored into the computational domain with a reverse sign following the method of image. 
Vertical force can be implemented in a straightforward way by replacing the Sinc function 
by its vertical derivative. The same interpolation function can be used for the extraction of 
the pressure wavefield at arbitrary receiver positions. The accuracy of the method of Hicks 
(2002) is illustrated in Figure 6a which shows a 3.75 Hz monochromatic wavefield computed 
in a homogeneous half space. The wave speed is 1500 m/s and the density is 1000 kg/m°. 
The grid interval is 100 m. The free surface is half a grid interval above the top of the FD 
grid and the method of image is used to implement the free surface boundary condition. The 
source is in the middle of the FD cell at 2 km depth. The receiver line is oriented in the Y 
direction. Receivers are in the middle of the FD cell in the horizontal plane and at a depth 
of 6 m just below the free surface. Comparison between the numerical and the analytical 
solutions at the receiver positions are first shown when the source is positioned at the closest 
grid point and the numerical solutions are extracted at the closest grid point (Figure 6b). The 
amplitude of the numerical solution is strongly overestimated because the numerical solution 
is extracted at a depth of 50 m below free surface (where the pressure vanishes) instead of 6 m. 
Second, a significant phase shift between numerical and analytical solutions results from the 
approximate positioning of the sources and receivers. In contrast, a good agreement between 
the numerical and analytical solutions both in terms of amplitude and phase is shown in 
Figure 6c where the source and receiver positioning is implemented with the windowed Sinc 
interpolation. 
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Fig. 6. a) Real part of a 3.75Hz monochromatic wavefield in a homogeneous half space. (b) 
Comparison between numerical (black) and analytical (gray) solutions at receiver positions 
when the closest grid point is used for both the source implementation and the extraction of 
the solution at the receiver positions on a coarse FD grid. (c) Same comparison between 
numerical (black) and analytical (gray) solutions at receiver positions when the Sinc 
interpolation with 4 coefficients is used for both the source implementation and the 
extraction of the solution at the receiver positions on a coarse FD grid. 


4. Realistic examples for acoustic and elastic propagations using FD formulations 


We shall provide two simple examples of seismic modeling using finite-differences methods 
both in the frequency and time approaches. The first example concerns seismic exploration 
problem where the acoustic approximation is often used while the second one is related 
to seismic risk mitigation where free surface effects including elastic propagation are quite 
important. 


4.1 3D EAGE/SEG salt model 


The salt model is a constant density acoustic model covering an area of 13.5 km x 13.5 km x 
4.2 km (Aminzadeh et al., 1997)(Figure 7). The salt model is representative of a Gulf Coast salt 
structure which contains salt sill, different faults, sand bodies and lenses. The salt model is 
discretized with 20 m cubic cells, representing an uniform mesh of 676 x 676 x 210 nodes. The 
minimum and maximum velocities in the Salt model are 1500 m/s and 4482 m/s respectively. 
We performed a simulation for a frequency of 7.33 Hz and for one source located at x = 
3600 m, y = 3.600 m and z = 100 m. The original model is resampled with a grid interval 
of 50 m corresponding to 4 grid points per minimum wavelength. The dimension of the 
resampled grid is 270 x 270 x 84 which represents 8.18 millions of unknowns after addition of 
the PML layers. Results of simulations performed with either in the frequency domain or in 
the time domain are compared in Figure 7. The time duration of the simulation is 15 s. 


We obtain a good agreement between the two solutions (Figure 7d) although we show a small 
phase shift between the two solutions at offsets greater than 5000 m. This phase shift results 
from the propagation in the high-velocity salt body. The direct-solver modeling is performed 
on 48 MPI process using 2 threads and 15 Gbytes of memory per MPI process. The memory 
and the elapsed time for the LU decomposition were 402 Gbytes and 2863 s, respectively. The 
elapsed time for the solution step for one right-hand side (RHS) is 1.4 s when we process 16 
RHS at a time during the solution step in MUMPS. The elapsed time for one time-domain 
simulation on 16 processors is 211 s. The frequency-domain approach is more than one order 
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Fig. 7. (a) Salt velocity model. (b-c) 7.33-Hz monochromatic wavefield (real part) computed 
with a finite-different formulation in the frequency domain (b) and in the time domain (c). 
(d) Direct comparison between frequency-domain (gray) and time-domain (black) solutions. 
The receiver line in the dip direction is: (top) at 100 m depth and at 3600 m depth in the cross 
direction. The amplitudes are corrected for 3D geometrical spreading. (bottom) at 2500 m 
depth and at 15000 m in the cross direction. 


of magnitude faster than the time-domain one when a large number of RHS members (2000) 
and a small number of processors (48) are used (Table 1). For a number of processors equal to 
the number of RHS members, the two approaches have the same cost. Of note, in the latter 
configuration (Np=N,1;), the cost of the two methods is almost equal in the case of the salt 
model (0.94 h versus 0.816 h). 


Over the last decades, simulations of wave propagation in complex media have been 
efficiently tackled with finite-difference methods (FDMs) and applied with success to 
numerous physical problems (Graves, 1996; Moczo et al., 2007). Nevertheless, FDMs suffer 
from some critical issues that are inherent to the underlying Cartesian grid, such as parasite 
diffractions in cases where the boundaries have a complex topography. To reduce these 
artefacts, the discretisation should be fine enough to reduce the ’stair-case’ effect at the 
free surface. For instance, a second-order rotated FDM requires up to 60 grid points per 
wavelength to compute an accurate seismic wavefield in elastic media with a complex 
topography (Bohlen & Saenger, 2006). Such constraints on the discretisation drastically 
restrict the possible field of realistic applications. Some interesting combinations of FDMs 
and finite-element methods (FEMs) might overcome these limitations (Galis et al., 2008). The 
idea is to use an unstructured FEM scheme to represent both the topography and the shallow 
part of the medium, and to adopt for the rest of the model a classical FDM regular grid. For the 
same reasons as the issues related to the topography, uniform grids are not suitable for highly 
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Model) Method |Pre. (hr) Sol. (hr) Total (hr)||Pre. (hr) Sol. (hr) Total (hr) 
Salt Time 0 39 39 0 0.94 0.94 
Salt |Frequency| 0.8 0.78 1.58 0.80 0.016 0.816 


Table 1. Comparison between time-domain and frequency-domain modeling for 32 (left) and 
2000 (right) processors. The number of sources is 2000. Pre. denotes the elapsed time for the 
source-independent task during seismic modeling (i.e., the LU factorization in the 
frequency-domain approach). Sol. denotes the elapsed time for multi-RHS solutions during 
seismic modeling (i.e., the substitutions in the frequency-domain approach). 
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Fig. 8. On the left, the French Riviera medium with complex topography and bathymetry: an 
hypothetical earthquake of magnitude 4.5 is at a depth of 10 km below the epicenter shown 
by a red ball. The simulation medium is 20 km by 20 km by 15 km. On the right, the related 
peak ground acceleration (PGA). Please note that the acceleration is always lower than one 
tenth of the Earth acceleration g 


heterogeneous media, since the grid size is determined by the shortest wavelength. Except in 
some circumstances, like mixing grids (Aoi & Fujiwara, 1999) or using non uniform Cartesian 
grids (Pitarka, 1999) in the case of a low velocity layer, it is almost impossible to locally 
adapt the grid size to the medium properties in the general case. From this point of view, 
FEMs are appealing, since they can use unstructured grids or meshes. Due to ever-increasing 
computational power, these kinds of methods have been the focus of a lot of interest and have 
been used intensively in seismology (Aagaard et al., 2001; Akcelik et al., 2003; Ichimura et al., 
2007). 


4.2 PGA estimation in the French Riviera 


Peak ground acceleration (PGA) are estimated using empirical attenuation laws calibrated 
through databases of seismic records of various areas: these laws should be adapted to each 
area around the world and European moderate earthquakes require a specific calibration 
(Berge-Thierry et al., 2003). Aside these attenuation laws, numerical tools as finite-differences 
time-domain methods allows the deterministic estimation of the peak ground acceleration 
(PGA) in specific areas of interest once the medium is known and the source specified. 
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Small areas as the French Riviera where a complex topography and bathymetric makes the 
simulation difficult. We would like to illustrate the procedure of time-domain simulation 
on this specific example (Cruz-Atienza et al., 2007). The Figure 8 shows a very simple model 
surrounding the city of Nice: the box is 20 km by 20 km by 15 km in depth. The P-wave velocity 
is 5700 m/s while the S-wave velocity is 3300 m/s and the density 2600 km/ m°. The water is 
characterized by a P-wave velocity of 1530 m/s while the density is about 1030 km/ m3. The 
grid step is 50 m and the time integration step is 0.005 s. 


The numerical simulation of an hypothetical earthquake of magnitude 4.5 at a depth of 10 km 
in the Mediterranean Sea provides us a deterministic estimation of the PGA as shown in the 
Figure 8. This small source is characterized upto a frequency of 3 Hz and we select a source 
time function with this expected spectral content. 


Successful applications have been proposed in the Los Angeles basin and is improved as we 
increase our knowledge about the medium of propagation and about the source location and 
its characterization. The PGA is estimated everywhere and one can see that increase of the 
PGA is observed at the sea bottom and nearby the coast. One can show that the amplification 
of PGA is decreased when considering the water layer at the expense of a longer duration of 
the seismic signal. 

Of course, various simulations should be performed using different models of the medium 
and for various source scenarii. These simulations could help to assess the variability of the 
acceleration for possible potential earthquakes and may be used for the mitigation of seismic 
risks. The importance of constraining the model structure should be emphasized and we can 
accumulate this knowledge through various and different initiatives performed for a more 
accurate reconstruction of the velocity structure (Rollet et al., 2002). One tool is the seismic 
imaging procedure we have underlined in this chapter. 


5. Finite-elements discontinuous Galerkin methods: a weak formulation 


Finite element methods, often more intensive in computer resources, introduce naturally 
boundary conditions in an explicit manner. Therefore, we expect improved accurate solutions 
with this numerical approach at the expense of computer requirements. The system of 
equations (14) in time has now a non-diagonal mass matrix while the system of equations 
(15) has a impedance matrix particularly ill-conditioned in 3D geometry taking into account 
its dimensionality. Therefore, for 2D geometries, the frequency formulation is still a quite 
feasible option while time domain approaches are there appealing when considering 3D 
geometries. Due to ever-increasing computational power, finite element methods using 
unstructured meshes have been the focus of increased interest and have been used extensively 
in seismology (Aagaard et al., 2001; Akcelik et al., 2003; Ichimura et al., 2007). 


Usually, the approximation order remains low, due to the prohibitive computational cost 
related to a non-diagonal mass matrix. However, this high computational cost can be avoided 
by mass lumping, a standard technique that replaces the large linear system by a diagonal 
matrix (Chin-Joe-Kong et al., 1999; Marfurt, 1984) and leads to an explicit time integration. 
Another class of FEMs that relies on the Gauss-Lobatto-Legendre quadrature points has 
removed these limitations, and allows for spectral convergence with high approximation 
orders. This high-order FEM, called the spectral element method (SEM) (Komatitsch & 
Vilotte, 1998; Seriani & Priolo, 1994) has been applied to large-scale geological models up 
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to the global scale (Chaljub et al., 2007; Komatitsch et al., 2008). The major limitation of 
SEM is the exclusive use of hexahedral meshes, which makes the design of an optimal mesh 
cumbersome in contrast to the flexibility offered by tetrahedral meshes. With tetrahedral 
meshes (Frey & George, 2008), it is possible to fit almost perfectly complex topographies 
or geological discontinuities and the mesh width can be adapted locally to the medium 
properties (h-adaptivity). The extension of the SEM to tetrahedral elements represents 
ongoing work, while some studies have been done in two dimensions on triangular meshes 
(Mercerat et al., 2006; Pasquetti & Rapetti, 2006). On the other hand, another kind of FEM 
has been proven to give accurate results on tetrahedral meshes: the Discontinuous Galerkin 
finite-element method (DG-FEM) in combination with the arbitrary high-order derivatives 
(ADER) time integration (Dumbser & Käser, 2006). Originally, DG-FEM has been developed 
for the neutron transport equation (Reed & Hill, 1973). It has been applied to a wide range 
of applications such as electromagnetics (Cockburn et al., 2004), aeroacoustics (Toulopoulos & 
Ekaterinaris, 2006) and plasma physics (Jacobs & Hesthaven, 2006), just to cite a few examples. 
This method relies on the exchange of numerical fluxes between adjacent elements. Contrary 
to classical FEMs, no continuity of the basis functions is imposed between elements and, 
therefore, the method supports discontinuities in the seismic wavefield as in the case of a 
fluid/solid interface. In such cases, the DG-FEM allows the same equation to be used for 
both the elastic and the acoustic media, and it does not require any explicit conditions on 
the interface (Kaser & Dumbser, 2008), which is, on the contrary, mandatory for continuous 
formulations, like the SEM (Chaljub et al., 2003). Moreover, the DG-FEM is completely 
local, which means that elements do not share their nodal values, contrary to conventional 
continuous FEM. Local operators make the method suitable for parallelisation and allow for 
the mixing of different approximation orders (p-adaptivity). 


5.1 3D finite-element discontinuous Galerkin method in the time domain 


Time domain approaches are quite attractive when considering explicit time integration. 
However, in most studies, the DG-FEM is generally used with high approximation orders. 
We present a low-order DG-FEM formulation with the convolutional perfectly matched layer 
(CPML) absorbing boundary condition (Komatitsch & Martin, 2007; Roden & Gedney, 2000) 
that is suitable for large-scale three-dimensional (3D) seismic wave simulations. In this 
context, the DG-FEM provides major benefits. 


The p-adaptivity is crucial for efficient simulations, in order to mitigate the effects of the very 
small elements that are generally encountered in refined tetrahedral meshes. Indeed, the 
p-adaptivity allows an optimised time stepping to be achieved, by adapting the approximation 
order according to the size of the elements and the properties of the medium. The benefit of 
such a numerical scheme is particularly important with strongly heterogeneous media. Due to 
the mathematical formulation we consider, the medium properties are assumed to be constant 
per element. Therefore, meshes have to be designed in such a way that this assumption 
is compatible with the expected accuracy. The discretization must be able to represent 
the geological structures fairly well, without over-sampling, while the spatial resolution of 
the imaging process puts constraints on the coarsest parameterisation of the medium. If 
we consider full waveform inversion (FWI) applications, the expected imaging resolution 
reaches half a wavelength, as shown by Sirgue & Pratt (2004). Therefore, following the 
Shannon theorem, a minimum number of four points per wavelength is required to obtain 
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such accuracy. These reasons have motivated the development of DG-FEM with low orders. 
We focus on the quadratic interpolation, which yields a good compromise between accuracy, 
discretisation and computational cost. 


5.1.1 3D time-domain elastodynamics 


It is worth to provide notations for specific manipulation of equations for DG-FEM 
approaches. The first-order hyperbolic system (8) under the so-called pseudo-conservative 
form can be written following the approach of Ben Jemaa et al. (2007) as 


p= Y dM) +f 


dE {x,y,z} 
Ade = }, (NoT) + Ad, (41) 
dE {x,y,z} 
with the definitions of the velocity and stress vectors as F = (vy Uy Vz)' and G = 


(01 02 03 Oxy Oxz Fyz)". Under this pseudo-conservative form, the RHS of (41) does not include 
any term that relates to the physical properties. The diagonal matrix A has been introduced 
in the system (8) and its inverse is required for the computation of the stress components 
(equation (41)). Matrices Mg and Nọ are constant real matrices (Etienne et al., 2010). The 
extension of the pseudo-conservative form for the visco-elastic cases could be considered with 
the inclusion of memory variables while the anisotropic case should be further analysed since 
the change of variable may depend on the physical parameters. Finally, in the equation (41), 
the medium density is denoted by p, while f and go are the external forces and the initial 
stresses, respectively. 


5.1.2 Spatial discretisation 


Following standard strategies of finite-element methods (Zienkiewicz et al., 2005), we want 
to approximate the solution of the equation (41) by means of polynomial basis functions 
defined in volume elements. The spatial discretisation is carried out with non-overlapping 
and conforming tetrahedra. We adopt the nodal form of the DG-FEM formulation (Hesthaven 
& Warburton, 2008), assuming that the stress and velocity vectors are approximated in the 
tetrahedral elements as follows 


dj 

T(x, t) = Bij (Xj, t) Pi (x) 
= 

= di 

Gj (X,t) = ) G(X), t) pi (x), (42) 
= 


where i is the index of the element, ¥ is the spatial coordinates inside the element, and t 
is the time. d; is the number of nodes or degrees of freedom (DOF) associated with the 
interpolating Lagrangian polynomial basis function gj relative to the j-th node located at 
position x je Vectors Ti j and fij are the velocity and stress vectors, respectively, evaluated at 
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Fig. 9. (a) Po element with one unique DOF. (b) P; element with four DOF. (c) P> element with 
10 DOF. 


the j-th node of the element. Although it is not an intrinsic limitation, we have adopted here 
the same set of basis functions for the interpolation of the velocity and the stress components. 
In the following, the notation P, refers to a spatial discretisation based on polynomial basis 
functions of degree k, and a P element is a tetrahedron in which a Pą scheme is applied. The 
number of DOF in a tetrahedral element is given by dj = (k+1)(k+2)(k+3)/6. For instance, 
in a Po element (Figure 9.a), there is only one DOF (the stress and velocity are constant per 
element), while in a P) element (Figure 9.b), there are four DOF located at the four vertices 
of the tetrahedron (the stress and velocity are linearly interpolated). It is worth noting that 
the Po scheme corresponds to the case of the finite-volume method (Ben Jemaa et al., 2009; 
2007; Brossier et al., 2008). For the quadratic approximation order P;, one node is added at the 
middle of each edge of the tetrahedron, leading to a total of 10 DOF per element (Figure 9.c). 
The first step in the finite-element formulation is to obtain the weak form of the elastodynamic 
system. To do so, we multiply the equation (41) by a test function ;, and integrate the system 
over the volume of the element i. For the test function, we adopt the same kind of function 
as used for the approximation of the solution. This case corresponds to the standard Galerkin 
method and can be written as 


| Pir POLE dV = | gir Y d(Moë)dV 
Vi Vi de {x,y,z} 


Lonnaëav = Por E (Wa Wren di), (43) 
Vi Vi de {x,y,z} 


where the volume of the tetrahedral element i is denoted by V;. For the purpose of clarity, 
we have omitted the external forces and stresses in the equation (43). Standard manipulations 
of finite-elements methods (integration by parts, Green theorem for fluxes along boundary 
surfaces) are performed as well as an evaluation of centered flux scheme for its non-dissipative 
property (Ben Jemaa et al., 2007; Delcourte et al., 2009; Remaki, 2000). Moreover, we assume 
constant physical properties per element. We define the tensorial product © as the Kronecker 


Modelling Seismic Wave Propagation for Geophysical Imaging 279 
product of two matrices A and B given by 
411B ... 41»B 
aeB=| . . . |, (44) 
aB i dumB 


where (n x m) denotes the dimensions of the matrix A. We obtain the expression 


s soal = 4 
pi(Z3 © Ki) = — V7 (Mo 8 Ei) + 5 L [(Px © Fix)Gj + (Pik 8 Gé] 
de{x,y,z} KEN; 
5 aat 5 a 
(A18 K;)d =- Yo (Wo8 Ep); + 5 2 [Qx @ Fix) + (Qik @ Gix)B| , (45) 
de{x,y,z} kEN; 


where Z3 represents the identity matrix. In the system (45), the vectors 0; and 0; should be 
red as the collection of all nodal values of the velocity and stress components in the element 
i. The system (45) indicates that the computations of the stress and velocity wavefields in one 
element require information from the directly neighbouring elements. This illustrates clearly 
the local nature of DG-FEM. The flux-related matrices P and Q are defined as follows 


P= Yo tice Mo 
de{x,y,z} 

Qx= YD rire No, 
de {x,y,z} 


where the component along the 0 axis of the unit vector 7i;, of the face S;, that points from 
element i to element k is denoted by n;,9, while we also introduce the mass matrix, the stiffness 
matrix and the flux matrices with 0 € {x,y,z} respectively, 


(Ki)j = i, Pir pj AV j,rE ll, di], 

(Eio) = for) gidV  j,rE [1 di], 

(Fik)rj = Í Pir PAS  j,r € [1, di] (46) 
Y Vik 

(Gix)rj = | Pir PAS rell di] jell dl. 
J Sik 


It is worth noting that, in the last equation of the system (46), the DOF of elements i and k 
appear (d; and d,, respectively) indicating that the approximation orders are totally decoupled 
from one element to another. Therefore, the DG-FEM allows for varying approximation 
orders in the numerical scheme. This feature is referred to as p-adaptivity. Moreover, given 
an approximation order, these matrices are unique for all elements (with a normalisation 
according to the volume or surface of the elements) and they can be computed before hand 
with appropriate integration quadrature rules. The memory requirement is therefore low, 
since only a collection of small matrices is needed according to the possible combinations of 
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approximation orders. The maximum size of these matrices is (dmax X dmax) Where dmax is 
the maximum number of DOF per element and the number of matrices to store is given by 
the square of the number of approximation orders mixed in the numerical domain. The four 
matrices Kj, Ei, Fj, and Gj, are computed by numerical integration using Hammer quadrature 
(Hammer & Stroud, 1958) and explicit forms of these matrices could be found in Etienne et al. 
(2010) for Po, P and P, orders. 


It should be mentioned that, in order to retrieve both the velocity and the stress components, 
the system (45) requires the computation of aah which can also be performed before hand. 
Note that, if we want to consider variations in the physical properties inside the elements, the 
pseudo-conservative form makes the computation of flux much easier and computationally 
more efficient than in the classical elastodynamic system. These properties come from the fact 
that, in the pseudo-conservative form, the physical properties are located in the left-hand side 
of the system (41). Therefore, no modification of the stiffness and flux matrices nor additional 
terms are needed in the system (45) to take into account the variation of properties. Only the 
mass matrix needs to be evaluated for each element and for each physical property according 
to the expression 


(Kiri = f xsl) girl) pA jr elt di, (47) 


where x;(%) represents the physical property (p; or one of the A; components) varying inside 
the element. 


5.1.3 Time discretisation 


The time integration of the system (45) relies on the second-order explicit leap-frog scheme 
that allows to compute alternatively the velocity and the stress components between a half 
time step. The system (45) can be written as 


n+} on—i 
‘ — D; 2 z 1 5 2, 
PT BK) = — E (MS) +5 L [(PxS Fat + (Pa & Gaag] 
be {x,y,z} kEN; 
go = g" n+l 
(A; BR = = D (No Q Eig) 3; 
ve 
+1 
+5 A [Qr Fi); na + (Qn D Ga)? ‘|, (48) 
2 KEN 


where the superscript n indicates the time step. We chose to apply the definition of the time 
step as given by Käser et al. (2008), which links the mesh width and time step as follows 
1 2r; 
At at 4 

<P ee Voi)” (49) 
where r; is the radius of the sphere inscribed in the element indexed by i, Vp; is the P-wave 
velocity in the element, and k; is the polynomial degree used in the element. Equation (49) is a 
heuristic stability criterion that usually works well. However, there is no mathematical proof 
for unstructured meshes that guarantees numerical stability. 
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Fig. 10. Speed-up observed when the number of MPI processes is increased from 1 to 256 for 
modelling with a mesh of 1.8 million P) elements. The ideal speed-up is plotted with a 
dashed line, the observed speed-up with a continuous line. These values were observed on a 


computing platform with bi-processor quad core Opteron 2.3 GHz CPUs interconnected with 
Infiniband at 20 Gb/s. 


5.1.4 Computational aspects 


The DG-FEM is a local method, and therefore it is naturally suitable for parallel computing. 
In our implementation, the parallelism relies on a domain-partitioning strategy, assigning 
one subdomain to one CPU. This corresponds to the single program mutiple data (SPMD) 
architecture, which means that there is only one program and each CPU uses the same 
executable to work on different parts of the 3D mesh. Communication between the 
subdomains is performed with the message passing interface (MPI) parallel environment 
(Aoyama & Nakano, 1999), which allows for applications to run on distributed memory 
machines. For efficient load balancing among the CPUs, the mesh is divided with the 
partitioner METIS (Karypis & Kumar, 1998), to balance the number of elements in the 
subdomains, and to minimise the number of adjacent elements between the subdomains. 
These two criteria are crucial for the efficiency of the parallelism on large-scale numerical 
simulations. Figure 10 shows the observed speed-up (i.e. the ratio between the computation 
time with one CPU, and the computation time with N CPUs) when the number of MPI 
processes is increased from 1 to 256, for strong scaling calculations on a fixed mesh of 
1.8 million P) elements. This figure shows good efficiency of the parallelism, of around 
80%. In our formulation, another key point is the time step, which is common for all of 
the subdomains. The time step should satisfy the stability condition given in equation (49) 
for every element. Consequently, the element with the smallest time step imposes its time 
step on all of the subdomains. We should mention here a more elaborate approach with 
local time stepping (Dumbser et al., 2007) that allows for elements to have their own time 
step independent of the others. Nevertheless, the p-adaptivity offered by DG-FEM allows 
mitigation of the computational burden resulting from the common time step as we shall see. 


5.1.5 Source excitation 


We proceed with the addition of the excitation to incremental increase of each involved field 
component. The excitation of a point source is projected onto the nodes of the element that 
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Fig. 11. (a) Cross-section of the mesh near the source position, indicated with a yellow star in 
the xy plane. This view represents the spatial support of the stress component in a Pp element 
containing the point source. (b) Same as (a) with a P; element. (c) Same as (a) with a P2 
element. 


0 


contains the source as follows 


an Gi (Xs) 
= s(t), (50) 
2A Pij (Xs) fy, Pi ()dV 


with s/ the nodal values vector associated to the excited component, t = nAt, ¥s the position 
of the point source and s(t) the source function. Equation (50) gives the source term that 
should be added to the right-hand side of equation (48) for the required components. It should 
be noticed that this term is only applied to the element containing the source. Depending 
on the approximation order, the spatial support of the source varies. Figure 11.a shows 
that the support of a Pp element is actually the whole volume of the element (represented 
on the cross-section with a homogeneous white area). In this case, no precise localisation 
of the source inside the element is possible due to the constant piece-wise interpolation 
approximation. On the other hand, in a P; element (Figure 11.b), the spatial support of the 
source is linear and allows for a rough localisation of the source. In a Py element (Figure 11.c), 
the quadratic spatial support tends to resemble the expected Dirac in space close to the source 
position. It should be noted that the limitations concerning source localisation also apply to 
the solution extraction at the receivers, according to the approximation order of the elements 
containing the receivers. 


5.1.6 Free surface condition 


Among the various approaches presented previously, we proceed by considering that the free 
surface follows the mesh elements. For the element faces located on the free surface, we use an 
explicit condition by changing the flux expression locally. This is carried out with the concept 
of virtual elements, which are exactly symmetric to the elements located on the free surface. 
Inside the virtual elements, we impose a velocity wavefield that is identical to the wavefield of 
the corresponding inner elements, and we impose an opposite stress wavefield on this virtual 
element. Thanks to the nodal formulation, the velocity is seen as continuous across the free 
surface, while the stress is equal to zero on the faces related to the free surface. 
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This is a quite natural approach similar to the one used in continuous finite-element methods 
where the test function is set to zero on the free surface boundary. 


5.1.7 Absorbing boundary condition 


We proceed through some simulations of wave propagation in a homogeneous, isotropic and 
purely elastic medium for an illustration of CPML conditions. The model size is 8 km x 8 
km x 8 km, and the medium properties are: Vp = 4000 m/s, Vs = 2310 m/s and p = 2000 
kg/ mê. An explosive source is placed at coordinates (x,= 2000 m, ys = 2000 m, zs = 4000 m) 
and a line of receivers is located at coordinates (3000 m < x; < 6000 m, y; = 2000 m, zr = 
4000 m) with 500 m between receivers. The conditions of the tests are particularly severe, 
since the source and the receivers are located close to the CPMLs (at a distance of 250 m), 
thus favouring grazing waves. The source signature is a Ricker wavelet with a dominant 
frequency of 3 Hz and a maximum frequency of about 7.5 Hz. Due to the explosive source, 
only P-wave is generated and the minimum wavelength is about 533 m. The mesh contains 
945,477 tedrahedra with an average edge of 175 m, making a discretisation of about 3 elements 
per Amin. Figures 12.c and 12.d show the results obtained with the P, interpolation and CPMLs 
of 10-elements width (L¢pmj = 1750 m) at all edges of the model. With the standard scale, no 
reflection can be seen from the CPMLs. When the amplitude is magnified by a factor of 100, 
some spurious reflections are visible. This observation is in agreement with the theoretical 
reflection coefficient (Reoeff = 0.1%) in equation (20). 

As shown by Collino & Tsogka (2001), the thickness of the absorbing layer plays an important 
role in the absorption efficiency. In Figures 12.a and 12.b, the same test was performed 
with CPMLs of 5-elements width (L¢pmj = 875 m) at all edges of the model. Compared to 
Figures 12.c and 12.d, the amplitude of the reflections have the same order of magnitude. 
Nevertheless, in the upper and left parts of the model, some areas with a strong amplitude 
appear close to the edges. These numerical instabilities arise at the outer edges of the CPMLs, 
and they expand over the complete model during the simulations. 


Instabilities of PML in long time simulations have been studied in electromagnetics 
(Abarbanel et al., 2002; Bécache et al., 2004). For elastodynamics, remedies have been 
proposed by Meza-Fajardo & Papageorgiou (2008) for an isotropic medium with standard 
PML. These authors proposed the application of an additional damping in the PML, onto 
the directions parallel to the layer, leading to a multiaxial PML (M-PML) which does not 
follow strictly the matching property of PML in the continuum and which has a less efficient 
absorption power. Through various numerical tests, Etienne et al. (2010) has shown that 
instabilities could be delayed outside the time window of simulation when considering 
extended M-PML from CPML. 


5.1.8 Saving computation time and memory 


Table 2 gives the computation times for updating the velocity and stress wavefields in 
one element for one time step, for different approximation orders, without or with the 
update of the CPML memory variables (i.e. elements located outside or inside the CPMLs). 
These computation times illustrate the significant increase with respect to the approximation 
order, and they allow an evaluation of the additional costs of the CPML memory variables 
computation from 40% to 60%. The effects of this additionnal cost have to be analysed in 
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Fig. 12. Snapshots at 1.6 s of the velocity component vy in the plane xy that contains the 
source location. CPMLs of 10-elements width are applied at all edges of the model. The 
modelling was carried out with P) interpolation. White lines, the limits of the CPMLs; black 
cross, the position of the source. (a) Real amplitude. (b) Amplitude magnified by a factor of 
100. (c) & (d) Same as (a) & (b) with CPMLs of 5-elements width. 


Approximation order Element outside CPML Element inside CPML 


Po 2.6 us 3.6 us 
Pi 5.0 us 8.3 us 
Py 21.1 us 29.9 us 


Table 2. Computation times for updating the velocity and stress wavefields in one element 
for one time step. These values correspond to average computation times for a computing 
platform with bi-processor quad core Opteron 2.3 GHz CPUs interconnected with Infiniband 
20 at Gb/s. 


the context of a domain-partitioning strategy. The mesh is divided into subdomains, using 
a partitioner. Figure 13.a shows the layout of the subdomains that were obtained with the 
partitioner METIS (Karypis & Kumar, 1998) along the xy plane used in the previous validation 
tests. The mesh was divided into 32 partitions, although only a few of these are visible on the 
cross-section in Figure 13.a. We used an unweighted partitioning, meaning that each partition 
contains approximately the same number of elements. 


The subdomains, partially located in the CPMLs, contain different numbers of CPML 
elements. In large simulations, some subdomains are totally located inside the CPMLs, and 
some others outside the CPMLs. In sucha case, the extra computation costs of the subdomains 
located in the absorbing layers penalise the whole simulation. Indeed, most of the subdomains 
spend 40% to 60% of the time just waiting for the subdomains located in the CPMLs to 
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Fig. 13. (a) Layout of the subdomains obtained with the partitioner METIS (Karypis & 
Kumar, 1998) along the xy plane that contains the source location. Grey lines, the limits of 
the CPMLs. The mesh was divided into 32 partitions, although only a few of these are visible 
on this cross-section. (b) View of the approximation order per element along the same plane. 
Black, the P elements; white, the P elements. 


complete the computations at each time step. For a better load balancing, we propose to 
benefit from the p-adaptivity of DG-FEM, using lower approximation orders in the CPMLs. 
Indeed, inside the absorbing layers, we do not need a specific accuracy, and consequently 
the approximation order can be decreased. Table 2 indicates that such a mixed numerical 
scheme is advantageous, since the computation time required for a Pp or P4 element located in 
the CPML is shorter than the computation time of a standard P; element. Figure 13.b shows 
the approximation order per element when P; is used in the CPMLs and P, in the rest of the 
medium. We should note here that the interface between these two areas is not strictly aligned 
to a cartesian axis, and has some irregularities due to the shape of the tetrahedra. Although it 
is possible to constrain the alignment of the element faces parallel to the CPML limits, we did 
not observe significant differences in the absorption efficiency whether the faces are aligned 
or not. 

Figure 14.a shows the seismograms computed when the modelling was carried out with Pz 
inside the medium and P, in the CPMLs. Absorbing layers of 10-elements width are applied 
at all edges of the model. For comparison, Figure 14.b shows the results obtained with Pz 
inside the medium and Pp in the CPMLs. In this case, the spurious reflections have significant 
amplitudes, preventing any use of these seismograms. On the other hand, the seismograms 
computed with the mixed scheme P2/P, show weak artefacts, and are reasonably comparable 
with the seismograms obtained with complete P) modelling. Therefore, taking into account 
that the computation time and the memory consumption of the P)/P, simulation are nearly 
half of those required with the full P) modelling, we can conclude that this mixed numerical 
scheme is of interest. It should be noticed that it is possible to adopt a weighted partitioning 
approach to overcome partly load balancing issues We should also stress that the saving in 
CPU time and memory provided with this kind of low-cost absorbing boundary condition is 
crucial for large 3D simulations, and this becomes a must in the context of 3D seismic imaging 
applications that require a lot of forward problems, such as FWI. 
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Fig. 14. (a) Seismograms of the velocity component vy. The amplitude of each seismogram is 
normalised. The modelling is done with P4 in the CPMLs and P; inside the medium. Black 
continuous line, numerical solution in large model without reflection in the time window; 
dashed line, numerical solution with 10-elements width CPMLs; grey line, residuals 
magnified by a factor of 10. (b) Same as (a) except the modelling is done with Pp in the 
CPMLs and P; inside the medium. 


5.1.9 Accuracy of DG-FEM with tetrahedral meshes 


There are a variety of studies in the literature concerning the dispersive and dissipative 
properties of DG-FEM with reference to wave-propagation problems. Let us quote few 
examples: Ainsworth et al. (2006) provided a theoretical study for the 1D case; Basabe et al. 
(2008) analysed the effects of basis functions on 2D periodic and regular quadrilateral meshes; 
and Kaser et al. (2008) discussed the convergence of the DG-FEM combined with ADER time 
integration and 3D tetrahedral meshes. More related to our particular concern here, Delcourte 
et al. (2009) provided a convergence analysis of the DG-FEM with a centred flux scheme and 
tetrahedral meshes for elastodynamics. They demonstrated the sensitivity of the DG-FEM to 
the mesh quality, and they proved that the convergence is limited by the second-order time 
integration we have used in the present study, despite the order of the basis function. Specific 
analysis of the convergence in the scheme we have presented could be found in Etienne et al. 
(2010). 


5.2 2D finite-element discontinuous Galerkin method in the frequency domain 


On land exploration seismology, there is a need to perform elastic wave modeling in area of 
complex topography such as foothills and thrust belts (Figure 15) in the frequency domain. 
Moreover, onshore targets often exhibit weathered layers with very low wave speeds in the 
near surface which require a locally-refined discretisation for accurate modeling. In shallow 
water environment, a mesh refinement is also often required near the sea floor for accurate 
modeling of guided and interface waves near the sea floor. Accurate modeling of acoustic and 
elastic waves in presence of complex boundaries of arbitrary shape and the local adaptation 
of the discretisation to local features such as weathered near surface layers or sea floor 
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Fig. 15. Application of the DG method in seismic exploration. (a) Velocity model 
representative of a foothill area affected by a hilly relief and a weathered layer in the near 
surface. (b) Close-up of the unstructured triangular mesh locally refined near the surface. (c) 
Example of monochromatic pressure wavefield. 


were two of our motivations behind the development of a discontinuous element method 
on unstructured meshes for acoustic and elastic wave modeling. 


5.2.1 hp-adaptive discontinuous Galerkin discretisation 


Similarly to the time formulation we adopt the nodal form of the DG formulation, assuming 
that the wavefield vector is approximated in triangular elements for 2D geometry which leads 
to the following expression, 


di 
da, x,y,z) = } tig (Ww, Xj, yjz) Pi l(w, xyz), (51) 
j=l 

where ü is the wavefield vector of components such as the following vector ii = (p, Vx, Vy, Vz) 
for acoustic propagation. The index of the element in an unstructured mesh is denoted by i. 
The expression ii;(w, x,y,z) denotes the wavefield vector in the element i and (x,y,z) are the 
coordinates inside the element i. In the framework of the nodal form of the DG method, 9;; 
denotes Lagrange polynomial and d; is the number of nodes in the element i. The position of 
the node j in the element i is denoted by the local coordinates (xj, yj, zj). 


In the frequency domain, the pseudo-conservative form (41) could be written in a 2D geometry 
as 
Mii = dg (Noti) +5, (52) 
de{x,y,z} 
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where Nii are linear fluxes and the source vector is denoted by §. Expressions of matrices M 
and N could be found in Brossier et al. (2010). 


The weak form of the system (52) is similar in the frequency domain and proceed by selecting 
a test function ;, and then an integration over the element volume V; which gives 


f Pir Mil; dV = f gir Y do (Noñ) dV + f PindidV, (53) 
Vi Vi de{x,y,z} Vi 


where the quantity r € [1,d;]. In the framework of Galerkin methods, we used the same 
function for the test function and the shape function. Similar procedures as for the 3D case 
and related to standard steps of the finite-element method lead to the discrete expression, 


(M;@K)#=- }, (o8 Eip) ii + ; L [(Qix Q Fix) di + (Qik © Gik) ir] + (Z@Kj) 3; 
de {x,y,z} KEN; 

(54) 
where the mass matrix X;, the stiffness matrix €; and the flux matrices F; and G; are similar to 
those defined for the 3D case (equation (46)). The matrix Q is also defined as for the 3D case 
(equation (46)) 


It is worth repeting that, in the equation (46), arbitrary polynomial order of the shape functions 
can be used in elements i and k indicating that the approximation orders are totally decoupled 
from one element to another. Therefore, the DG allows for varying approximation orders in 
the numerical scheme, leading to the p-adaptivity. 


The equation (54) can be recast in matrix form as 


Bu=s. (55) 


5.2.2 Which interpolation orders to choose? 


For the shape and test functions, we used low-order Lagrangian polynomials of orders 0, 1 
and 2, referred to as Pk, k € 0,1,2 in the following (Brossier, 2009; Etienne et al., 2009). Let 
us remind that our motivation behind seismic modeling is to perform seismic imaging of the 
subsurface by full waveform inversion, the spatial resolution of which is half the propagated 
wavelength and that the physical properties of the medium are piecewise constant per 
element in our implementation of the DG method. The spatial resolution of the FWI and the 
piecewise constant representation of the medium direct us towards low-interpolation orders 
to achieve the best compromise between computational efficiency, solution accuracy and 
suitable discretisation of the computational domain. The Po interpolation (or finite volume 
scheme) was shown to provide sufficiently-accurate solution on 2D equilateral triangular 
mesh when ten cells per minimum propagated wavelength are used (Brossier et al., 2008), 
while 10 cells and 3 cells per propagated wavelengths provide sufficiently-accurate solutions 
on unstructured triangular meshes with the P; and the P) interpolation orders, respectively 
(Brossier, 2011). Of note, the Py scheme is not convergent on unstructured meshes when 
centered fluxes are used (Brossier et al., 2008). This prevents the use of the Py scheme in 
3D medium where uniform tetrahedral meshes do not exist (Etienne et al., 2010). A second 
remark is that the finite volume scheme on square cells is equivalent to second-order accurate 
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FD2D DG? DG?P DGP 
nq 1 1 3 6 
nl 9 5-9 13-25 24-48 


Table 3. Number of nodes per element (ng) and number of non-zero coefficients per row of 
the impedance matrix (nz) for the FD and DG methods. The number nz depends on the 
number of wavefield components involved in the r.h.s of the first-order wave equation nger- 


FD stencil (Brossier et al., 2008) which is consistent with a discretisation criterion of 10 grid 
points per wavelength (Virieux, 1986). Use of interpolation orders greater than 2 would allow 
us to use coarser meshes for the same accuracy but these coarser meshes would lead to an 
undersampling of the subsurface model during imaging. On the other hand, use of high 
interpolation orders on mesh built using a criterion of 4 cells par wavelength would provide 
an unnecessary accuracy level for seismic imaging at the expense of the computational cost 
resulting from the dramatic increase of the number of unknowns in the equation (55). 


The computational cost of the LU decomposition depends on the numerical bandwidth of 
the matrix, the dimension of the matrix (i.e., the number of rows/columns) and the number of 
non-zero coefficients per row (nz). The dimension of the matrix depends in turn of the number 
of cell (ne), of the number of nodes per cell (ng) and the number of wavefield components 
(Nwave) (ranging from 3 to 5 in 2D geometry). The number of nodes in a 2D triangular element 
is given by Hesthaven & Warburton (2008) and leads to the following expression ng = (k + 
1)(k+ 2) /2 where k denotes the interpolation order similar to what is done in the 3D geometry. 


The numerical bandwidth is not significantly impacted by the interpolation order. The 
dimension of the matrix and the number of non-zero elements per row of the impedance 
matrix are respectively given by fwave X ng X Mce and (1 + Mneigh) X Na X Nder + 1, where 
Nneigh is the number of neighbor cell (3 in 2D geometry) and nger is the number of wavefield 
components involved in the rh.s of the velocity-pressure wave equation, equation (52). 
Table 3 outlines the number of non-zero coefficients per row for the mixed-grid FD and DG 
methods. Increasing the interpolation order will lead to an increase of the number of non-zero 
coefficients per row, a decrease of the number of cells in the mesh and an increase of the 
number of nodes in each element. The combined impact of the 3 parameters nz, ho, Ng 
on the computational cost of the DG method makes difficult the definition of the optimal 
discretisation of the frequency-domain DG method. The medium properties should rather 
drive us towards the choice of a suitable discretisation. 


One must underline that the LU factorization is quite demanding in computer memory and 
has also some drawbacks for scalability, suggesting that nodes with high memory should be 
preferred at the expense of the CPU numbers. 


5.2.3 Boundary conditions and source implementation 


Absorbing boundary conditions are implemented with unsplitted PML in the 
frequency-domain DG method (Brossier, 2011) following the same approach than for 
the FD method: one can see that the PML implementation in the frequency is straightforward. 
We have found that constraining the meshing to have edges of elements in the PML zone 
parallel to the direction of dissipation of the waves improves the efficiency. 
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Free surface boundary condition is implemented with the method of image. A virtual 
cell is considered above the free surface with the same velocity and the opposite pressure 
components to those below the free surface. This allows us to fulfill the zero pressure 
condition at the free surface while keeping the correct numerical estimation of the particle 
velocity at the free surface. Using these particle velocities and pressures in the virtual cell, 
the pressure flux across the free surface interface vanishes, while the velocity flux is twice 
the value that would have been obtained by neglecting the flux contribution above the free 
surface. As in the FD method, this boundary condition has been implemented by modifying 
the impedance matrix accordingly without introducing explicitely the virtual element in the 
mesh. The rigid boundary condition is implemented following the same principle except that 
the same pressure and the opposite velocity are considered in the virtual cell. 

Concerning the source excitation, the point source at arbitrary positions in the mesh is 
implemented by means of the Lagrange interpolation polynomials for k > 1. This means 
that the source excitation is performed at the nodes of the cell containing the source with 
appropriate weights corresponding to the projection of the physical position of the source on 
the polynomial basis. When the source is located in the close vicinity of a node of a triangular 
cell, all the weights are almost zero except that located near the source. In the case of the P2 
interpolation, a source close to the vertex of the triangular cell is problematic because the 
integral of the P basis function over the volume of the cell is zero for nodes located at the 
vertex of the triangle. In this case, no source excitation will be performed (see equation (54)). 
To overcome this problem specific to the P, interpolation, one can use locally a P; interpolation 
in the element containing the source at the expense of the accuracy or distribute the source 
excitation over several elements or express the solution in the form of local polynomials (i.e., 
the so-called modal form) rather than through nodes and interpolating Lagrange polynomials 
(i.e., the so-called nodal form). 

Another issue is the implementation of the source in Po equilateral mesh. If the source 
is excited only within the element containing the source, a checker-board pattern is 
superimposed on the wavefield solution. This pattern results from the fact that one cell out of 
two is excited in the DG formulation because the DG stencil does not embed a staggered-grid 
structure (the unexcited grid is not stored in staggered-grid FD methods; see Hustedt et al. 
(2004) for an illustration). To overcome this problem, the source can be distributed over 
several elements of the mesh or P; interpolation can be used in the area containing the sources 
and the receivers, while keeping Po interpolation in the other parts of the model (Brossier 
et al., 2010). 

Of note, use of unstructured meshes together with the source excitation at the different nodes 
of the element contribute to mitigate the checker-board pattern in the in P; and P) schemes. 
The same procedure as for the source is used to extract the wavefield solution at arbitrary 
receiver positions. 


6. Realistic examples for highly contrasted and strongly heterogeneous media 
using finite-elements methods 


We shall consider two examples for the illustration of the Discontinuous Galerkin approach. 
The first one is related to the problem of 3D wave propagation inside an active volcano 
using the time-domain approach while the second one deals with the problem of 2D wave 
propagation above a oil reservoir using the frequency-domain approach. 
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6.1 The volcano La Soufriére 
6.1.1 Characteristics of the model 


La Soufriére of Guadeloupe (France) is one of nine active volcanoes of Lesser Antilles. It 
belongs to a recent volcanic system situated in the south part of the Basse-Terre. A P-wave 
velocity model of the volcano has been obtained by first arrival time tomography (Coutant 
et al., 2010). Figure 16 is the reconstructed Vp velocity model that reveals the existence 
of a high velocity zone below the dome of La Soufrière. The dimensions of the model are 
1400 m x 1400 m x 1000 m in xyz respectively. We consider a constant Poisson ratio of 0.25 
to assess the S-wave velocity model from Vp. The velocity ranges from 660 m/s to 3800 m/s 
for Vp and 380 m/s to 2200 m/s for Vs. Considering a maximum frequency of 25 Hz, the 
minimum wavelength is about 15 m. In addition, we consider a constant density equal to 
2000 kg/m. Absorbing layers of CPML type with a thickness of 300 m are added at each 
side edge of the model as well at the bottom edge. Therefore, the complete dimensions of the 
numerical model are 2000 m x 2000 m x 1300 m. 


Fig. 16. Topography of the volcano La Soufriére with the underlying reconstructed Vp 
velocity structure. The position of the dynamite shot is indicated with a yellow circle and the 
receivers with black triangles. 


6.1.2 Construction of the tetraedral mesh 


The mesh has been built with the mesher TETGEN (Si, 2006) combined with an iterative 
h-refinement procedure to obtain a locally adapted mesh to the velocity field(with an average 
of 3 elements per minimum wavelength Anin): a cross-section is shown in the left panel of 
the Figure 17. For building this mesh, we have started our iterative reconstruction with a 
uniform mesh shown in the right panel of the Figure 17. After the sixth refinement iteration, 
the discretization criteria are met. Areas of high velocities are correlated with the parts of the 
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Modelling time 5s 
Nb elements 4.6 million 
Nb unknowns 414 million 
Min/Max element edge |1.29 - 58.62 m 
Nb time steps 37 787 
Nb CPUs 512 
Total memory 10 GB 
Memory per CPU 14-23 MB 
Computation time 6h45 min. 
Time / unknown / time step} 0.79 us 


Table 4. Statistics of the modelling for La Soufriére performed on an IBM Blue Gene machine. 
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Fig. 17. On the left, cross-section of the P-wave velocity model in the plane xz in the middle 
of the volcano La Soufriére. The back line represents the topography. On the right, initial 
mesh of the volcano La Soufriére from which we deduce automatically the one used for 
modeling by adapting the mesh size to the local P-wave velocity. Absorbing layers of CPML 
type with a thickness of 300 m are added at each side edge of the model. 
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mesh where the elements are the largest ones. On the contrary, near the free surface, we find 
the finest elements. 


6.1.3 Numerical result 


We have performed 3D simulations with the Discontinuous Galerkin Finite-Element Method 
in the time domain. The computations have been performed on a Blue Gene machine with 
512 processors. The statistics for these computations are given in Table 4. 


The configuration of the seismic acquisition is given in Figure 16. This is a quasi-2D system 
with a profile according to the East-West direction, which includes 100 single-component 
receivers (vz) with 10 m between receivers. The source is a shot of dynamite. For the numerical 
simulations, we used an explosive source with a Ricker function of dominant frequency of 10 
Hz (maximum frequency 25 Hz). We present in Figure 18 a comparison between the observed 
and computed data. Despite significant uncertainties and approximations (source function, 
Poisson ratio, density, absence of attenuation, low signal to noise ratio), there are striking 
similarities in the data. In particular, the seismic traces exhibit well marked discontinuities 
related to the strong velocity contrasts and the complex topography of the volcano La Soufrière. 
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Fig. 18. (a) Recorded seismograms (component vz). (b) Computed seismograms. Some 
similarities between both set of data are highlighted with color shapes. 
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Fig. 19. The synthetic Valhall model for (a) P-wave and (b) S-wave velocities. Panel (c) 
represents a zoom of the shallow mesh. 


6.2 Application 2D in the frequency-domain: the synthetic Valhall application 


This 2D application is based on a synthetic representation of the Valhall zone in the North Sea, 
Norway. This model is representative of oil and gas fields in shallow water environments of 
the North Sea (Munns, 1985). The model is described as an heterogeneous P- and S- wave 
velocity model (Figure 19a-b). The water layer is only 70 m depth. The main targets are 
a gas cloud in the large sediment layer, and in a deeper part of the model, the trapped oil 
underneath the cap rock, which is formed of chalk. Gas clouds are easily identified by the low 
P-wave velocities, whereas their signature is much weaker in the Vs model, as gas does not 
affect S-waves propagation. 
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Fig. 20. Frequency-domain data for the hydrophone component at 4 Hz. The data (real-part) 
are plotted in the source/receiver domain. 


In order to investigate seismic imaging in such environment, the selected acquisition mimics 
a four-component ocean-bottom cable survey (Kommedal et al., 2004), as is deployed on the 
field. A line of 315 explosive sources is positioned 5 m below water surface to simulated 
air-gun sources and a cable of 315 3-components sensors is located on the sea floor (1 
hydrophone and 2 geophones). This geological setting, which is composed of a significant 
soft sea-bed with high Poisson’ratio due to soft and unconsolidated sediments leads to a 
particularly ill-posed problem for S-wave velocity reconstruction, due to the relatively small 
shear-wave velocity contrast at the sea bed, which prevents recording of significant P-to-S 
converted waves. 


For the meshing of the model, the narrow velocity range in most parts of the model requires 
the use of a regular mesh as much as possible for computational efficiency. However, to 
correctly discretize the shallow-water layer and liquid-solid interface, a p-adapted mesh 
implemented with a mixed Po-P, interpolation is chosen (Figure 19c): a refined unstructured 
P, layer of cells is used for the first 130 m of the subsurface for accurate modeling of the 
interface waves at the liquid/solid interface, and for accurate positioning of the sources, 
located 5 m below the surface, and of the receivers, located on the sea floor. Below this 130 m 
depth zone, a regular equilateral mesh is used in combination of a Po interpolation. 


Figure 20 illustrates an example of frequency-domain data (real-part of the complex-value 
wavefield) at 4Hz. These data are plotted in the source/receiver domain for the full 
acquisition survey. The diagonal part of the figure represents the collocation of the source 
and the closest receiver, as the source moves in the acquisition. The Figure 21 illustrates 
a time-domain shot-gather for the 3 components of the sensors and a source located at 
position 4 km : the frequency-domain solutions at all the receiver positions have been 
computed for the single source at all the frequencies of the source spectrum, between 0 and 


Modelling Seismic Wave Propagation for Geophysical Imaging 295 


b) Offset (km) C) 


Time (s) 
Time (s) 
Time (s) 


Vx component Vz component 


Fig. 21. Time-domain shot-gather for a source located at position 4 km. The data are 
computed from frequency-domain simulation and Fourier transformed in the time-domain 
for (a) the horizontal geophone, (b) the vertical geophone and (c) the hydrophone 
components. 


13 Hz. These frequency-domain complex-values data have then been Fourier transform to the 
time-domain. The time-domain show specific properties of propagation in such environments 
: the hydrophone and the vertical geophones are mainly sensitive to P-wave arrivals that 
dominate the elastic propagation in soft-sediment zones. The horizontal geophone allows 
however to record some late P-to-S conversions, which could be used to image the Vs model 
from seismic imaging methods. 


7. Conclusion 


We have presented mainly two families of techniques for solving partial differential equations 
for elastodynamics: some finite-differences formulations in both time domain and frequency 
domain and some finite-element methods also in both time domain and frequency domain. 
Both approaches have appealing features, especially when considering seismic imaging where 
numerous forward problems should be performed. Such classification helps to understand 
the advantages and limitations of each particular method to model a specific physical 
phenomenon 


The discretization of the strong formulation of the partial differential equations has been 
presented through finite-difference techniques. These approaches are easy to implement 
and quite flexible. They are currently the methods of choice for large-scale modelling 
and inversion in exploration geophysics, especially in the marine environment. They may 
however demand a very fine discretization when the earth model contains large contrasts; 
and accurately modelling the responses around a sharp interface is quite challenging. We 
have introduced various perspectives as summation-by-parts formulation or the immersed 
boundary approach as well as simple mesh deformation might broaden the use of 
finite-differences techniques by avoiding the stair-case approximation. 
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The weak formulation expressed in the finite-element methods has been considered under 
the specific family of Discontinuous Galerkin approaches. The use of test functions gives us 
more freedom and the integral form provides us flexibility in the meshing. However, they 
lead to numerical challenges: they are more difficult to implement than the finite-difference 
method, they are often more expensive in computational time and memory, and they are 
more complicated to use because the accuracy of the response depends on the quality of the 
meshing. Therefore, they are not intensively used for seismic imaging and are until now more 
oriented to seismic modeling in the final reconstructed model. 


It should be noticed that attempts exist to combine the advantages of these methods in one 
approach for computing elastic fields, at least for specific applications. Even, one can think 
that decoupling the inverse problem procedure and the forward problem is possible: we 
can flip-flop between the two forward problem formulations inside iterations of the inverse 
problem. 


When the modelling method serves as the kernel of an inversion algorithm, additional 
constrains generally appear because the gradient of the misfit functional needs to be 
evaluated. The choice of the modelling approach notably depends (1) on the needed accuracy, 
(2) the efficiency in evaluating the solution and the gradient of the misfit functional in an 
inversion algorithm, and (3) the simplicity of use. 


Finally, the practical implementation shall probably be adapted to the data acquisition. 
Densely sampled acquisition in exploration geophysics with or without blending, or in 
lithospheric investigation with the recent deployment of sensors such as the USarray 
experiment challenges our modelling choice. This seems to indicate that development in 
modelling and associated inversion approaches remain crucial to improve our knowledge 
of the subsurface, notably by extracting more information from the, ever larger, recorded data 
sets. 
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